1 Preprocessing data

sceset is the SingleCellExperiment object that contains all the single cell data

sceset <- readRDS("../../scFT-paper_rds/20191016data_integration_all_sceset.rds")
dim(sceset)
## [1] 22110  4330
## Add secretory annotations
secretory <- readRDS("../../man_analysis2_20180629/rds/20190120Fresh_secretory_9clusters_clincluster.rds")

secretory$cell_type[secretory$clincluster_final == "C10"] <- "Cell cycle (C9)"
secretory$cell_type[secretory$clincluster_final == "C3"] <- "Differentiated (C3)"
secretory$cell_type[secretory$clincluster_final == "C4"] <- "KRT17 (C4)"
secretory$cell_type[secretory$clincluster_final == "C6"] <- "Stress (C6)"
secretory$cell_type[secretory$clincluster_final == "C8"] <- "EMT (C7)"
secretory$cell_type[secretory$clincluster_final == "C9"] <- "Immune (C8)"
secretory$cell_type[secretory$clincluster_fina %in% c("C1","C2","C5")] <- "Queiscent"
secretory$secretory_subtype <- secretory$cell_type
cancer_secretory <- secretory ## the clustering results of cancer_secretory
rm(secretory)

1.1 Removing doublets

sceset$secretory_subtype <- cancer_secretory$cell_type[match(colnames(sceset), colnames(cancer_secretory))]

seurat <- as.Seurat(sceset)
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')
seurat <- NormalizeData(seurat)
seurat <- ScaleData(object = seurat)
seurat <- FindVariableFeatures(seurat, selection.method = "vst", nfeatures = 2000)
seurat <- RunPCA(seurat)
seurat <- RunUMAP(seurat, dims = 1:15, seed.use = 12345)
## Warning: Cannot add objects with duplicate keys (offending key: UMAP_),
## setting key to 'umap_'
seurat <- RunTSNE(seurat, dims = 1:15, seed.use = 12345)

seurat <- FindNeighbors(seurat, dims = 1:15)
seurat <- FindClusters(seurat, resolution = 0.6)
# DimPlot(seurat, reduction = "tsne", label = T)

seurat$disease <- NA
seurat$disease[seurat$Patient2 %in% c(11543, 11545, 11553)] <- "ovarian"
seurat$disease[seurat$Patient2 %in% c(15066, 15091, 15072)] <- "endometrial"
seurat$disease[seurat$Patient2 %in% c(33572, 33778, 34350, 34659,35773)] <- "benign"

sweep.res.list_kidney <- paramSweep_v3(seurat, PCs = 1:15, sct = FALSE)
sweep.stats_kidney <- summarizeSweep(sweep.res.list_kidney, GT = FALSE)
bcmvn_kidney <- find.pK(sweep.stats_kidney)

## Homotypic Doublet Proportion Estimate 
annotations <- seurat$seurat_clusters ## Clustering results
homotypic.prop <- modelHomotypic(annotations)       
nExp_poi <- round(0.01*length(seurat@active.ident))  ## 1% doublet rate
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))

## Run DoubletFinder with varying classification stringencies
seurat <- doubletFinder_v3(seurat, PCs = 1:10, pN = 0.25, pK = 0.005, nExp = nExp_poi, reuse.pANN = FALSE)
seurat <- doubletFinder_v3(seurat, PCs = 1:10, pN = 0.25, pK = 0.005, nExp = nExp_poi.adj, reuse.pANN = "pANN_0.25_0.005_43")

seurat@meta.data$DF_hi.lo <- seurat@meta.data$DF.classifications_0.25_0.005_43
seurat@meta.data$DF_hi.lo[which(seurat@meta.data$DF_hi.lo == "Doublet" & seurat@meta.data$DF.classifications_0.25_0.005_36 == "Singlet")] <- "Doublet_lo"
seurat@meta.data$DF_hi.lo[which(seurat@meta.data$DF_hi.lo == "Doublet")] <- "Doublet_hi"

# DimPlot(seurat, reduction = "tsne", group.by="DF_hi.lo", plot.order=c("Doublet_hi","Doublet_lo","Singlet"), colors.use=c("black","gold","red"))

colData(sceset)$DF_hi.lo <- seurat$DF_hi.lo
# all(colnames(seurat)==colnames(sceset))

benign_ft <- seurat[,seurat$DF_hi.lo == "Singlet" & seurat$type=="Benign"]
# saveRDS(benign_ft, "rds/20191030benignFT_1857cells_SeuratObj.rds")

In total, 1857 cells from benign FT were left for analysis.

dim(benign_ft)
## [1] 22110  1857

1.2 UMAP plots

# p1 <-  DimPlot(benign_ft, reduction = "TSNE", group.by = "Patient2", pt = 0.4) + 
#     ggtitle(paste("t-SNE plot \n Benign FT (n = ", ncol(benign_ft), ")", sep = ""))
DimPlot(benign_ft, reduction = "UMAP", group.by = "Patient2", pt = 0.4) + 
    ggtitle(paste("UMAP plot \n Benign FT (n = ", ncol(benign_ft), ")", sep = ""))

# cowplot::plot_grid(p1,p2)
# ggsave("../revision_plots/SI/FigS3_benignFT_UMAPtSNE_colByPatient.png", width = 10, height = 4)
FeaturePlot(benign_ft, reduction = "UMAP",features = c("PTPRC","COL1A1","EPCAM","KRT7","PAX8","CAPS"), ncol  = 3, pt.size = 0.4)

# ggsave("../revision_plots/SI/FigS3_benignFT_UMAPcolourByMarkers.png", width = 14, height = 8)

2 Filter out non-secretory cells

From the cells from benign donors we filtered out non-secretory cells.

# sum(benign_ft$seurat_clusters %in% c(8,10:12) | benign_ft$FACS == "CD45+EPCAM-") # 122 leukocytes
secretory <- benign_ft[,!benign_ft$seurat_clusters %in% c(8,10:12) & benign_ft$FACS != "CD45+EPCAM-"] # remove CD45+ cells

# sum(secretory$FACS == "STROMAL") # 11 stroma cells
secretory <- secretory[, !secretory$FACS %in% c("STROMAL")]

dim(secretory)
## [1] 22110  1724

1724 FTE cells in total

# sum(secretory$seurat_clusters[secretory$type== "Benign"] %in% c(3, 7)) # 515 ciliated cells
secretory <- secretory[, !secretory$seurat_clusters %in% c(3, 7)] # remove cilicated cells

cpm <- log2(calculateCPM(secretory@assays$RNA@counts) + 1)
secretory <- secretory[, (cpm["KRT7",] > 0.5 | cpm["PAX8",] > 0.5) &
                          cpm["EPCAM",] > 0.5  &
                          cpm["PTPRC",] == 0 
                       ]
rm(cpm)
table(secretory$Patient2)
## 
## 33572 33778 34350 34659 35773 
##   210    53   130   356    16
# 33572 33778 34350 34659 35773 
#   210    53   130   356    16 
secretory_new <- secretory # secretory from the benign donors
# saveRDS(secretory, "rds/20191030benign_765secretory_seuratObj.rds")
cancer_secretory <- as.Seurat(cancer_secretory)
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')
cancer_secretory$type <- "Cancer"
secretory <- merge(secretory_new, cancer_secretory)

secretory$secretory_subtype[is.na(secretory$secretory_subtype) & secretory$type == "Cancer"] <- "Unclassified"

# saveRDS(secretory, "rds/seurat_secretory_merged20191007.rds")

3 Batch effects

Before we integrated the secretory cells from cancer and benign patients, let’s first see the batch effects between them.

secretory <- NormalizeData(secretory)
secretory <- ScaleData(object = secretory)
secretory <- FindVariableFeatures(secretory, selection.method = "vst", nfeatures = 2000)
secretory <- RunPCA(secretory)
ElbowPlot(secretory, 40)

secretory <- RunUMAP(secretory, dims = 1:14, seed.use = 123345)
secretory <- RunTSNE(secretory, dims = 1:14, seed.use = 123345)

secretory <- FindNeighbors(secretory, dims = 1:14, seed.use = 123345)
secretory <- FindClusters(secretory, resolution = 0.6)

It is clear that batch effects exist between the two batches as shown in the figure below. So the next step we did is to remove the batch effects.

p1 <- DimPlot(secretory, reduction = "umap", label = T)
## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.
p2 <- DimPlot(secretory, reduction = "umap", group.by = "type")
cowplot::plot_grid(p1,p2)

3.1 Data integration (standard workflow)

This section integrated different batches.

seu_list <- list()

seu_list[[1]] <- secretory[,secretory$type == "Cancer"]
seu_list[[1]] <- seu_list[[1]][,seu_list[[1]]$secretory_subtype != "Immune (C8)" & seu_list[[1]]$secretory_subtype != "Unclassified"]

seu_list[[2]] <- secretory[,secretory$Patient2 %in% c("34350","33778","35773","33572")]
seu_list[[3]] <- secretory[,secretory$Patient2 == "34659"]

for(itor in 1:3){
    seu_list[[itor]] <- NormalizeData(seu_list[[itor]], verbose = FALSE)
    seu_list[[itor]] <- FindVariableFeatures(seu_list[[itor]], 
                                             selection.method = "vst", 
                                             nfeatures = 1000, verbose = FALSE)
}

n_dim <- 15
ftsc.anchors <- FindIntegrationAnchors(object.list = seu_list, 
                                       dims = 1:n_dim,
                                       k.anchor = 5) 
## Computing 2000 integration features
## Scaling features for provided objects
## Finding all pairwise anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 1779 anchors
## Filtering anchors
##  Retained 1167 anchors
## Extracting within-dataset neighbors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 1645 anchors
## Filtering anchors
##  Retained 1085 anchors
## Extracting within-dataset neighbors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 1454 anchors
## Filtering anchors
##  Retained 1156 anchors
## Extracting within-dataset neighbors
ftsc.integrated <- IntegrateData(anchorset = ftsc.anchors, dims = 1:n_dim)
## Merging dataset 3 into 2
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data
## Merging dataset 1 into 2 3
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data
DefaultAssay(ftsc.integrated) <- "integrated"

# Run the standard workflow for visualization and clustering
ftsc.integrated <- ScaleData(ftsc.integrated, verbose = FALSE)
ftsc.integrated <- RunPCA(ftsc.integrated, npcs = n_dim, verbose = FALSE)
ftsc.integrated <- RunUMAP(ftsc.integrated, reduction = "pca", dims = 1:n_dim, seed.use = 12345)
p1 <- DimPlot(ftsc.integrated, reduction = "umap", group.by = "type", pt.size = 0.5)
p2 <- DimPlot(ftsc.integrated, 
              reduction = "umap", 
              group.by = "secretory_subtype", pt.size = 0.5)
cowplot::plot_grid(p1,p2, rel_widths = c(1,1.1))

p1 <- DimPlot(ftsc.integrated[,ftsc.integrated$type == "Cancer" &
                                  !is.na(ftsc.integrated$secretory_subtype)], 
              reduction = "umap", group.by = "secretory_subtype") + ggtitle("Cancer")
p2 <- DimPlot(ftsc.integrated[,ftsc.integrated$type == "Benign"],
              group.by = "type", reduction = "umap") + ggtitle("Benign")
cowplot::plot_grid(p1,p2)

3.2 Clustering results

ftsc.integrated <- FindNeighbors(ftsc.integrated, dims = 1:n_dim, seed.use = 12345, k.param = 5)
## Computing nearest neighbor graph
## Computing SNN
ftsc.integrated <- FindClusters(ftsc.integrated, resolution = 0.6)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1453
## Number of edges: 18343
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7525
## Number of communities: 7
## Elapsed time: 0 seconds
# ftsc.integrated <- readRDS("../../scFT-paper_rds/ftsc.integrated_20191013.rds")

DimPlot(ftsc.integrated, label = T, reduction = "umap")

table(ftsc.integrated$secretory_subtype[ftsc.integrated$type == "Cancer"], ftsc.integrated$seurat_clusters[ftsc.integrated$type == "Cancer"])
##                      
##                         0   1   2   3   4   5   6
##   Cell cycle (C9)      10   2   1   3   0   0   7
##   Differentiated (C3)  21   1  65  80   0   0   0
##   EMT (C7)              3   2   0   0   1  34   0
##   KRT17 (C4)           14  10   1  20 139   2   0
##   Stress (C6)         161  22  72  11   4   0   2
ftsc.integrated$secretory_subtype_benign <- NA
ftsc.integrated$secretory_subtype_benign[ftsc.integrated$type == "Benign" & ftsc.integrated$seurat_clusters %in% c(4)] <- "KRT17"

ftsc.integrated$secretory_subtype_benign[ftsc.integrated$type == "Benign" & ftsc.integrated$seurat_clusters %in% c(9)] <- "Cell cycle"

ftsc.integrated$secretory_subtype_benign[ftsc.integrated$type == "Benign" & ftsc.integrated$seurat_clusters %in% c(5)] <- "EMT"

ftsc.integrated$secretory_subtype_benign[ftsc.integrated$type == "Benign" & ftsc.integrated$seurat_clusters %in% c(2,3)] <- "Differentiated"


ftsc.integrated$secretory_subtype_benign[ftsc.integrated$type == "Cancer" & ftsc.integrated$secretory_subtype == "KRT17 (C4)"] <- "KRT17"

ftsc.integrated$secretory_subtype_benign[ftsc.integrated$type == "Cancer" & ftsc.integrated$secretory_subtype == "Cell cycle (C9)"] <- "Cell cycle"

ftsc.integrated$secretory_subtype_benign[ftsc.integrated$type == "Cancer" & ftsc.integrated$secretory_subtype == "EMT (C7)"] <- "EMT"

ftsc.integrated$secretory_subtype_benign[ftsc.integrated$type == "Cancer" & ftsc.integrated$secretory_subtype == "Differentiated (C3)"] <- "Differentiated"

3.3 Scoring the cells from benign samples

new_markers <- readRDS("../../scFT-paper_rds/20190213new_markers.rds")
ftsc.integrated$KRT17_score <- colSums(ftsc.integrated@assays$RNA@data[new_markers$C4,])
ftsc.integrated$CC_score <- colSums(ftsc.integrated@assays$RNA@data[new_markers$C10,])
ftsc.integrated$EMT_score <- colSums(ftsc.integrated@assays$RNA@data[new_markers$EMT,])
ftsc.integrated$Diff_score <- colSums(ftsc.integrated@assays$RNA@data[new_markers$C3,])
ftsc.sceset <- as.SingleCellExperiment(ftsc.integrated)

df.plot <- cbind(ftsc.sceset@reducedDims$UMAP, 
                 ftsc.sceset@colData[,c("KRT17_score", "CC_score", "EMT_score", "Diff_score", "type", "secretory_subtype_benign")])
df.plot <- as.data.frame(df.plot)
# df.plot <- df.plot[!is.na(df.plot$secretory_subtype_benign),]
df.plot1 <- df.plot[df.plot$type == "Benign",]
df.plot2 <- df.plot[df.plot$type == "Cancer",]

scaleScore <- function(x, max.limit = 3) {
    x <- scale(x, center = T, scale = T)
    x[x > max.limit] <- max.limit
    return(x)
}

df.plot1$KRT17_score <- scaleScore(df.plot1$KRT17_score)
df.plot1$CC_score <- scaleScore(df.plot1$CC_score)
df.plot1$Diff_score <- scaleScore(df.plot1$Diff_score)
df.plot1$EMT_score <- scaleScore(df.plot1$EMT_score)

df.plot2$KRT17_score <- scaleScore(df.plot2$KRT17_score)
df.plot2$CC_score <- scaleScore(df.plot2$CC_score)
df.plot2$Diff_score <- scaleScore(df.plot2$Diff_score)
df.plot2$EMT_score <- scaleScore(df.plot2$EMT_score)


plotScoreUmap <- function(dat, column){
        ggplot(data=dat, aes_string(x="UMAP_1", y="UMAP_2", fill = column)) + 
        geom_point(alpha = 0.5, pch = 21, col = "grey80") +
        theme_classic() +
        scale_fill_gradient2(midpoint=0,
                              low="grey90", mid="white",high="red", space ="Lab" ) +
        theme(plot.title = element_text(size=14, face="bold", hjust = 0.5))
}

plist <- list()
plist[[1]] <- plotScoreUmap(dat = df.plot1, column = "KRT17_score") + ggtitle("Benign (KRT17 score)") 
plist[[2]] <- plotScoreUmap(dat = df.plot1, column = "CC_score") + ggtitle("Benign (cell cycle score)")
plist[[3]] <- plotScoreUmap(dat = df.plot1, column = "Diff_score") + ggtitle("Benign (differentiated score)")
plist[[4]] <- plotScoreUmap(dat = df.plot1, column = "EMT_score") + ggtitle("Benign (EMT score)")
plist[[5]] <- plotScoreUmap(dat = df.plot2, column = "KRT17_score") + ggtitle("Cancer (KRT17 score)")
plist[[6]] <- plotScoreUmap(dat = df.plot2, column = "CC_score") + ggtitle("Cancer (cell cycle score)")
plist[[7]] <- plotScoreUmap(dat = df.plot2, column = "Diff_score") + ggtitle("Cancer (differentiated score)")
plist[[8]] <- plotScoreUmap(dat = df.plot2, column = "EMT_score") + ggtitle("Cancer (EMT score)")

cowplot::plot_grid(plotlist = plist, ncol = 4)

# ggsave("../revision_plots/SI/benign_cancer_subtype_score_plots.png", height = 8, width = 18)
# ggsave("../revision_plots/SI/benign_cancer_subtype_score_plots.tiff", height = 8, width = 18)

# ggsave("../revision_plots/SI/benign_cancer_subtype_score_plots.pdf", height = 8, width = 18)

The expression of marker genes is expressed in the genes at the similar location.

# saveRDS(ftsc.integrated, "rds/ftsc.integrated_20191013.rds")

4 Technical

sessionInfo()
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.4
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] fields_9.6                  maps_3.3.0                 
##  [3] spam_2.2-2                  dotCall64_1.0-0            
##  [5] modes_0.7.0                 ROCR_1.0-7                 
##  [7] gplots_3.0.1.1              KernSmooth_2.23-15         
##  [9] DoubletFinder_2.0.1         reshape2_1.4.3             
## [11] scales_1.0.0                Seurat_3.0.2               
## [13] dplyr_0.7.8                 edgeR_3.24.3               
## [15] limma_3.38.3                scater_1.10.1              
## [17] ggplot2_3.2.0.9000          SingleCellExperiment_1.4.1 
## [19] SummarizedExperiment_1.12.0 DelayedArray_0.8.0         
## [21] BiocParallel_1.16.5         matrixStats_0.54.0         
## [23] Biobase_2.42.0              GenomicRanges_1.34.0       
## [25] GenomeInfoDb_1.18.1         IRanges_2.16.0             
## [27] S4Vectors_0.20.1            BiocGenerics_0.28.0        
## 
## loaded via a namespace (and not attached):
##  [1] Rtsne_0.15               ggbeeswarm_0.6.0        
##  [3] colorspace_1.4-0         ggridges_0.5.1          
##  [5] XVector_0.22.0           listenv_0.7.0           
##  [7] npsurv_0.4-0             ggrepel_0.8.0           
##  [9] codetools_0.2-16         splines_3.5.2           
## [11] R.methodsS3_1.7.1        lsei_1.2-0              
## [13] knitr_1.21               jsonlite_1.6            
## [15] ica_1.0-2                cluster_2.0.7-1         
## [17] png_0.1-7                R.oo_1.22.0             
## [19] sctransform_0.2.0        HDF5Array_1.10.1        
## [21] compiler_3.5.2           httr_1.4.0              
## [23] assertthat_0.2.0         Matrix_1.2-15           
## [25] lazyeval_0.2.1           htmltools_0.3.6         
## [27] tools_3.5.2              rsvd_1.0.2              
## [29] bindrcpp_0.2.2           igraph_1.2.3            
## [31] gtable_0.2.0             glue_1.3.0              
## [33] GenomeInfoDbData_1.2.0   RANN_2.6.1              
## [35] Rcpp_1.0.0               gdata_2.18.0            
## [37] ape_5.2                  nlme_3.1-137            
## [39] DelayedMatrixStats_1.4.0 gbRd_0.4-11             
## [41] lmtest_0.9-36            xfun_0.4                
## [43] stringr_1.4.0            globals_0.12.4          
## [45] irlba_2.3.3              gtools_3.8.1            
## [47] future_1.14.0            zlibbioc_1.28.0         
## [49] MASS_7.3-51.1            zoo_1.8-4               
## [51] rhdf5_2.26.2             RColorBrewer_1.1-2      
## [53] yaml_2.2.0               reticulate_1.10         
## [55] pbapply_1.4-0            gridExtra_2.3           
## [57] stringi_1.2.4            caTools_1.17.1.1        
## [59] bibtex_0.4.2             Rdpack_0.10-1           
## [61] SDMTools_1.1-221         rlang_0.4.0             
## [63] pkgconfig_2.0.2          bitops_1.0-6            
## [65] evaluate_0.13            lattice_0.20-38         
## [67] purrr_0.3.0              Rhdf5lib_1.4.2          
## [69] bindr_0.1.1              labeling_0.3            
## [71] htmlwidgets_1.3          cowplot_0.9.4           
## [73] tidyselect_0.2.5         plyr_1.8.4              
## [75] magrittr_1.5             R6_2.3.0                
## [77] pillar_1.3.1             withr_2.1.2             
## [79] fitdistrplus_1.0-14      survival_2.43-3         
## [81] RCurl_1.95-4.11          tsne_0.1-3              
## [83] tibble_2.0.1             future.apply_1.3.0      
## [85] crayon_1.3.4             plotly_4.8.0            
## [87] rmarkdown_1.11           viridis_0.5.1           
## [89] locfit_1.5-9.1           data.table_1.12.0       
## [91] metap_1.1                digest_0.6.18           
## [93] tidyr_0.8.2              R.utils_2.7.0           
## [95] munsell_0.5.0            beeswarm_0.2.3          
## [97] viridisLite_0.3.0        vipor_0.4.5
# R version 3.5.2 (2018-12-20)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS Mojave 10.14.4
# 
# Matrix products: default
# BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
# [1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#  [1] class_7.3-15                DoubletFinder_2.0.1         reshape2_1.4.3              scales_1.0.0               
#  [5] Seurat_3.0.2                dplyr_0.7.8                 edgeR_3.24.3                limma_3.38.3               
#  [9] scater_1.10.1               ggplot2_3.2.0.9000          SingleCellExperiment_1.4.1  SummarizedExperiment_1.12.0
# [13] DelayedArray_0.8.0          BiocParallel_1.16.5         matrixStats_0.54.0          Biobase_2.42.0             
# [17] GenomicRanges_1.34.0        GenomeInfoDb_1.18.1         IRanges_2.16.0              S4Vectors_0.20.1           
# [21] BiocGenerics_0.28.0        
# 
# loaded via a namespace (and not attached):
#   [1] backports_1.1.3          plyr_1.8.4               igraph_1.2.3             lazyeval_0.2.1          
#   [5] splines_3.5.2            listenv_0.7.0            usethis_1.4.0            digest_0.6.18           
#   [9] htmltools_0.3.6          viridis_0.5.1            gdata_2.18.0             magrittr_1.5            
#  [13] memoise_1.1.0            cluster_2.0.7-1          ROCR_1.0-7               remotes_2.0.2           
#  [17] globals_0.12.4           RcppParallel_4.4.3       R.utils_2.7.0            prettyunits_1.0.2       
#  [21] colorspace_1.4-0         ggrepel_0.8.0            xfun_0.4                 callr_3.1.1             
#  [25] crayon_1.3.4             RCurl_1.95-4.11          jsonlite_1.6             bindr_0.1.1             
#  [29] survival_2.43-3          zoo_1.8-4                ape_5.2                  glue_1.3.0              
#  [33] gtable_0.2.0             zlibbioc_1.28.0          XVector_0.22.0           leiden_0.3.1            
#  [37] pkgbuild_1.0.2           Rhdf5lib_1.4.2           future.apply_1.3.0       HDF5Array_1.10.1        
#  [41] bibtex_0.4.2             Rcpp_1.0.0               metap_1.1                viridisLite_0.3.0       
#  [45] reticulate_1.10          rsvd_1.0.2               SDMTools_1.1-221         tsne_0.1-3              
#  [49] htmlwidgets_1.3          httr_1.4.0               gplots_3.0.1.1           RColorBrewer_1.1-2      
#  [53] ica_1.0-2                pkgconfig_2.0.2          R.methodsS3_1.7.1        uwot_0.1.3              
#  [57] locfit_1.5-9.1           tidyselect_0.2.5         labeling_0.3             rlang_0.4.0             
#  [61] munsell_0.5.0            tools_3.5.2              cli_1.0.1                devtools_2.0.1          
#  [65] ggridges_0.5.1           stringr_1.4.0            yaml_2.2.0               npsurv_0.4-0            
#  [69] fs_1.2.6                 processx_3.2.1           knitr_1.21               fitdistrplus_1.0-14     
#  [73] caTools_1.17.1.1         purrr_0.3.0              RANN_2.6.1               bindrcpp_0.2.2          
#  [77] pbapply_1.4-0            future_1.14.0            nlme_3.1-137             R.oo_1.22.0             
#  [81] compiler_3.5.2           rstudioapi_0.9.0         curl_3.3                 beeswarm_0.2.3          
#  [85] plotly_4.8.0             png_0.1-7                lsei_1.2-0               tibble_2.0.1            
#  [89] stringi_1.2.4            ps_1.3.0                 desc_1.2.0               RSpectra_0.13-1         
#  [93] lattice_0.20-38          Matrix_1.2-15            pillar_1.3.1             Rdpack_0.10-1           
#  [97] lmtest_0.9-36            RcppAnnoy_0.0.11         data.table_1.12.0        cowplot_0.9.4           
# [101] bitops_1.0-6             irlba_2.3.3              gbRd_0.4-11              R6_2.3.0                
# [105] KernSmooth_2.23-15       gridExtra_2.3            vipor_0.4.5              sessioninfo_1.1.1       
# [109] codetools_0.2-16         MASS_7.3-51.1            gtools_3.8.1             assertthat_0.2.0        
# [113] pkgload_1.0.2            rhdf5_2.26.2             rprojroot_1.3-2          withr_2.1.2             
# [117] sctransform_0.2.0        GenomeInfoDbData_1.2.0   grid_3.5.2               tidyr_0.8.2             
# [121] DelayedMatrixStats_1.4.0 Rtsne_0.15               ggbeeswarm_